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Abstract Severe population declines can reduce species to small populations, offering permissive conditions for 
deleterious processes. For example, following such events, species can become prone to inbreeding and genetic drift 
which can lead to a loss of genetic diversity and evolutionary potentials. Hynobius chinensis is a poorly studied very 
rare and declining endangered amphibian species endemic to China in Changyang County. We investigated adult 
census population size by monitoring breeding populations from 2015 to 2018, developed microsatellite markers 
from the transcriptome and used them to investigate genetic diversity, and a population bottleneck in this species. 
We found H. chinensis in 4 different localities in a total area of 2.18 km’ and estimated the overall adult census 
population size at 386—404 individuals. The adult census size (mean + SE) per breeding pond ranged from 44 + 6 
to 141 + 8 individuals and appeared smaller than that reported in closely related species in undisturbed habitats. We 
developed and characterized 13 microsatellite markers in total. Analysis of data at 7 loci (N = 118) in Hardy-Weinberg 
equilibrium gathered from the largest population showed that genetic diversity level was low. The average number 
of alleles per locus was 2.14. The observed and expected heterozygosities averaged 0.38 and 0.40, respectively. The 
inbreeding coefficient was —0.06. All tests performed to investigate a population bottleneck, i.e. The Garza-Williamson 
test, Heterozygosity excess test, Mode shift test of allele frequency, and effective population size estimates detected 
a population bottleneck. The contemporary and the historical effective population sizes were estimated at 36 and 234 
individuals, respectively. We argue that as bottleneck effects, the studied population may have become prone to genetic 
drift and inbreeding, losing microsatellite alleles and heterozygosity. Our results suggest that populations of H. chinensis 
may have been extirpated in the study area. 
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It can also help formulate appropriate conservation 
strategies (Amos and Balmford, 2001; Frankham, 2005). 
Theoretical and empirical works show that severe 


1. Introduction 


Elucidating factors and processes affecting genetic 


diversity in species of conservation concern can help — population declines or population bottlenecks can affect 


understand the ecology and evolution of such species. 
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genetic diversity in several ways (Amos and Balmford, 
2001; Frankham, 2005; Hoffmann et al., 2017; Willi and 
Hoffmann, 2009). For examples, by reducing species 
to small populations, population bottlenecks can reduce 
the standing genetic diversity and impair the ability 
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to fix novel adaptive genetic mutations (Frankham, 
2005; Hoffmann et al., 2017). Still, these demographic 
events can offer permissive conditions for inbreeding 
and random drift, which can cause a continuous loss of 
genetic diversity and of individual fitness or inbreeding 
depression (Amos and Balmford; Frankham, 2005; 
Heredia-Bobadilla et a/., 2017; Jordan et al., 2009; Spear 
et al., 2006; Storfer et al., 2014; Willi and Hoffmann, 
2009). Due to these effects, affected species can lose 
their ability to respond to novel selection pressures and 
adapt novel environmental conditions (Frankham, 2005; 
Hoffmann et al., 2017; Willi and Hoffmann, 2009). 
Population bottlenecks can also render species vulnerable 
to extinction due to stochastic and catastrophic events 
(Amos and Balmford, 2001; Frankham, 2005; Willi and 
Hoffmann, 2009). 

Severe population bottlenecks cause the values of 
different population genetics parameters, including 
heterozygosity and allele frequency at selectively neutral 
loci such as microsatellites to differ from those expected 
at Hardy-Weinberg equilibrium across several generations 
(Garza and Williamson, 2001; Luikart et al., 1998; Peery 
et al., 2012). Different tests based on the analysis of 
multilocus microsatellite data within populations allow 
detecting severe population bottlenecks that occurred 
in the recent history of species by detecting violation of 
theoretical expectations under this equilibrium (Garza 
and Williamson, 2001; Luikart et al., 1998; Peery et al., 
2012). These include the Mode shift test, the Garza- 
Williamson test (the updated version of the Garza 
M-ratio test), and tests of Heterozygosity excess (Garza 
and Williamson, 2001; Luikart et al., 1998; Peery et al., 
2012). 

Potential factors causing population bottlenecks 
and loss of genetic diversity may vary between 
organisms (Rivera-Ortíz et al., 2015). For amphibians, 
anthropogenic destruction and degradation of natural 
habitats is probably among the most important (Allentoft 
and O’Brien, 2010; Rivera-Ortiz et al., 2015). Studies 
that analyze microsatellite data in amphibian species 
whose natural habitats have been considerably reduced 
and degraded by human activities show that such species 
have undergone fragmentation in gene pools. They have 
reduced census and effective sizes and exhibit molecular 
signs of recent population bottlenecks, inbreeding and 
loss of microsatellite alleles and heterozygosity (Allentoft 
and O’Brien, 2010; Heredia-Bobadilla et al., 2017; 
McCartney-Melstad and Shaffer, 2015; Rivera-Ortíz et al., 
2015; Storfer et al., 2014; Sugawara et al., 2015). 

The Chinese salamander, Hynobius chinensis (Caudata, 


Hynobiidae) is an amphibian species endemic to China 
(Fei and Ye, 2016; Frost 2018). Classified as Endangered 
by the IUCN Red List (IUCN, 2018), literature shows 
that H. chinensis is threatened by anthropogenic habitat 
degradation and destruction, causing habitat loss for 
populations (Fei and Ye, 2016). As effects of these threats, 
this species has undergone population declines and is 
very rare (Fei and Ye, 2016; IUCN, 2018). It is restricted 
to mountain forests in Changyang County in Yichang 
(Hubei province) -Note that Changyang County covers 
3 412 km’ (Source: Wikipedia)- at altitudes ranging from 
1 400 to 1 500 m (Fei and Ye, 2016; Frost, 2018). 

H. chinensis is a poorly known and studied species. 
Data on many aspects of Biology, including landscape 
use, genetic diversity, current population size, and 
demographic history are scarce. H. chinensis has been 
mostly involved in studies addressing systematic and 
taxonomic issues (Fei and Ye, 2016; Frost, 2018; Wang 
et al., 2007). 

The nature of threats facing H. chinensis, the level 
of rarity, and the declining population trend reported 
in the literature let suspecting that this species might 
have experienced severe population bottlenecks and 
lost genetic diversity in the recent past (Apodaca et al., 
2012; Heredia-Bobadilla et al., 2017; Sugawara et al., 
2016; Sugawara et al., 2015; Storfer et al., 2014; Wang 
et al., 2017). If the history of H. chinensis conforms to 
this hypothesized trajectory, we should expect to see 
representative populations of this species exhibiting 
molecular signs of population bottlenecks detectable by 
microsatellite markers, low levels of genetic diversity and 
reduced sizes (Apodaca et al., 2012; Heredia-Bobadilla 
et al., 2017; Storfer et al., 2014; Sugawara et al., 2016; 
Sugawara et al., 2015; Wang et al., 2017). 

We investigated (1) adult census population size, 
developed microsatellite markers and used them to 
investigate (ii) genetic diversity, and (iii) molecular 
signatures of a population bottleneck in H. chinensis. 
Our aim was to contribute to our understanding of the 
demography and genetic viability of this endangered 
species. 


2. Materials and Methods 


2.1. Study area and field procedure Published literature 
shows that H. chinensis has one breeding period per year 
which begins in November (Fei and Ye, 2016). During 
the breeding period, adult individuals aggregate, mate 
and reproduce in ponds (Fei and Ye, 2016). An individual 
female lays eggs in a pair of sacs, with an egg sac's length 
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measuring tens of cm (Fei and Ye, 2016). 

From November 2015 to April 2016, we searched 
for H. chinensis individuals and egg sacs in ponds at 
altitudes ranging from 1300 to 1600 m in Changyang (Fei 
and Ye, 2016; Frost, 2018). Egg sacs were considered 
as a precursory sign of the presence of H. chinensis as 
this species is the only Hynobiid salamander distributed 
in Changyang (Fei and Ye, 2016). We found breeding 
populations of H. chinensis in at different sites of different 
localities. 

Data from H. nebulosus showed that the number of 
adult male individuals could be accurately determined 
by examining breeding aggregations because a male 
individual would spend several weeks in a breeding pond 
(Kusano, 1980). We examined a breeding population of 
H. chinensis at a breeding pond found in one locality: 
Dongtoutang. The pond was small and shallow (Table 
1). And, the pond's water was enough clear to observe 
objects present in the pond. Based on our previous 
field experience, and on information provided by local 
people, we were expecting to observe the largest breeding 
population at this pond. At the landscape scale, the area 
was at 5096 covered by a degraded natural deciduous 
forest of Quercus sp. (Fagaceae). The pond's surrounding 
area was covered by various herbs (mainly Poaceae, 
Oxalidaceae, and Polygonaceae). The minimum distance 
to a forest patch was 25 m. 

H. chinensis individuals were carefully collected 
from the pond, using a net (Heyer et al., 1994). The 
salamanders were sexed based on the morphology of the 
cloacae (Fei and Ye, 2016), marked by toe-clipping and 
released (Heyer et al., 1994). The number of captured 
individuals was recorded for each sex. After a time 
period ranging from of 1 to 4 weeks, the population was 
re-examined. For each sex, the number of “captured”, 
i.e. non-toe-clipped, and the number of “recaptured” 
i.e. toe-clipped individuals were recorded (Heyer et а/., 
1994). Pairs of egg sacs were progressively individually 
tagged as they were spawned and their total number was 
determined at the end of the breeding season (Chen et al., 
2016; Kusano, 1980). The investigation was finished on 
March 25, 2016. There was no additional novel pair of 
egg sacs since March 2; and no unmarked male individual 
since March 8, 2016. Then, the number of adult females, 
F and the number of adult males, M were estimated as the 
number of pairs of egg sacs (Chen et al., 2016; Kusano 
and Miyashita, 1984; Kusano, 1980) and toe-clipped male 
individuals, respectively (Heyer ef al., 1994). The sex 
ratio, SR was estimated as SR=M/F and the adult census 
population size, №, as N, =M+F. 
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To minimize the impact on the populations 
(McCartney-Melstad and Shaffer, 2015), as the value of 
SR was close to that reported in closely related species 
(Kusano, 2012), for the breeding populations found in 
other localities and for the remaining study period, the 
values of N, were indirectly estimated. In summary, for 
each population, the number of adult females, F, was 
estimated based on egg clutch number as previously 
explained. And, N,, was estimated as: №, = F, + (F,*SR), 
where л represents a given population. For all the 
localities, the values of №, was determined for 2015—2016; 
2016—2017 and 2017—2018 breeding seasons, at the same 
ponds and following this method. 

From 22 different egg clutches, 22 H. chinensis 
hatched larvae samples were randomly collected. All these 
larvae and the toe clips were kept in anhydrous (98%) 
ethanol at —20 °С in the laboratory for DNA-based 
experiments. 

During the fieldwork, to identify threats facing H. 
chinensis in the study area, local people were interviewed 
for a better understanding of the evolution of landscapes 
and of the usage of the species by humans. The species? 
mating behavior was occasionally observed at the 
breeding pond found in Dongtoutang. 


2.2. Development of microsatellite markers An adult 
H. chinensis individual was brought to the laboratory, 
euthanized by injection of MS-22 (3g/l), and dissected. 
A tissue sample was taken from the liver. Total RNA was 
extracted (Omega Bio-Tek, USA) with TRIzol Regent® 
(Invitrogen, CA, USA). The purity, concentration and 
integrity of RNA were assessed by Nanodrop, Qubit 2.0 
and Agilent 2100 bioanalyzer, respectively. Using lug of 
total RNA as input, a library for Illumina sequencing was 
prepared following the NEBNext Ultra™ RNA protocol. 
The library preparation was sequenced on an Illumina 
HiSeq 2500 platform (pair-end reads, 125 bp). 

Sequencing results were analyzed and reads containing 
adaptors, reads containing poly-N, and low-quality reads 
(i.e. with more than 50% О < 10 bases) were filtered out. 
And, the clean reads (i.e. the remaining reads) were De 
novo assembled, using Trinity v2.0.6 software (Grabherr 
et al., 2011). 

The Simple Sequence Repeat Identification Tool 
available at http://archive. gramene.org/db/markers/ 
ssrtool/ was used to identify microsatellite loci in the 
assembled H. chinensis’ unigene DNA sequences. 
The following searching criteria were applied: (i) the 
minimum number of repeats: 5 times; (11) the repeat type: 
tetranucleotide; and (111) flanking sequence length: = 18 
bp. Using the software Primer v3.0, 74 microsatellite 
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loci were screened and 119 primer pairs were designed, 
applying the following criteria: (1) primer length 18—23 
bp, with 20 bp as optimum; (2) product size (bp): 100— 
250 bp; (3) primer melting temperature: 55-69 °C , with 
60 “С as optimum and (4) GC% content: 40-70%, with 
5096 as optimum, and (5) GC clump: 2. 

Genomic DNA was extracted from all the toe tissue 
and larvae samples, using TIANamp Genomic Kit 
(TIANGEN). The 119 primer pairs were synthesized 
(Applied Biosystems, Beijing, China) and tested by 
gradient PCR on the standard agarose gel. Individual 
PCR amplification was performed as following: an initial 
denaturation at 94 °С for 5 min was followed by 35 
cycles of denaturation at 94 °С and at a primer's specific 
annealing temperature for 30 sec, and an extension at 72 
°С for 45 sec and finally at 72 °С for 8 min. PCR was 
performed in a reaction volume of 10 ul containing 0.6 pl 
of template DNA, 5 ul of Taq DNA polymerase (Premix), 
of 0.3 ul of forward primer, 0.3 pl of reverse primer, and 
3.8 ul of ddH,O. In total, 51 primer pairs consistently 
amplified at optimal annealing temperature and produced 
clear bands in the expected size range. These primer pairs 
were retained for polymorphism test. 

The 5' end of each forward primer was labeled with 
one of the fluorescent dyes FAM, TAMRA or HEX, and 
polymorphism was tested, using data from the 22 larvae 
samples. In summary, PCR amplification of the 51 primer 
pairs was performed, following the protocol above 
described. For each individual larva, genotyping at the 51 
microsatellite loci was performed with an ABI Prism 3730 
Genetic Analyzer (Applied Biosystems, Beijing, China). 
Microsatellite alleles were scored, using GeneScan 
ROX 400 size standard and GeneMapper v4.0 software 
(Applied Biosystems, Beijing, China). The existence of 
genotyping errors (i.e. stuttering, large allele dropout) and 
of null alleles within the microsatellite data (N — 22) was 
checked, using the software Micro-Checker v2.2.3 (Van 
Oosterhout et al., 2004). 

An exact test for Hardy-Weinberg equilibrium 
(HWE) for each locus and a test for potential linkage 
disequilibrium (LD) between all pairs of loci were 
implemented in the online Genepop software (Raymond 
and Rousset, 1995), using Bonferroni-corrected P-values 
(Rice, 1989). Exact P-values were determined by 
Markov Chain method, using the following settings: 
dememorization: 10 000; batches: 100; iteration per batch: 
1 000 (Raymond and Rousset, 1995). The values of basic 
genetic diversity indices, exclusion probabilities, and 
polymorphism information content (PIC) were determined 
for each locus and for overall loci, using Cervus v3.0 


software (Kalinowski et al., 2007). Genotypes data were 
analyzed and the sampled H. chinensis adult individuals 
were genotyped at polymorphic loci in HWE, without LD, 
null alleles and/or genotyping errors (Table 2), following 
the procedure above described. 


2.3. Analysis of genetic diversity and bottleneck effects 
For the adult samples’ genotype data, the observed 
heterozygosity (H,), expected heterozygosity (Н), and 
the number of alleles (4) were calculated for each locus 
and for the overall loci, using Arlequin v3.5 software 
(Excoffier and Lischer, 2010). The inbreeding coefficient, 
F, was calculated for each locus and for the overall loci, 
using FSTAT v2.9.3 software (Goudet, 2001). 

Different approaches were used to investigate 
genetic signatures of a population bottleneck. First, the 
signature of a population bottleneck was tested by Garza- 
Williamson test. In summary, Garza-Williamson index, 
GW was estimated for each locus as GW=k/(r+1), where 
k is the number of alleles, and r the allele size range 
(Excoffier and Lischer, 2010). And, the overall GW 
value was compared to the theoretical critical value for 
rejecting the null hypothesis of demographic stability 
that is 0.68 (Garza and Williamson, 2001; Peery et al., 
2012). The value of GW is expected to be below this 
critical value in bottlenecked populations. The core idea 
behind Garza-Williamson test is that during population 
bottleneck, reduction in k is proportionally important 
than the reduction in r. Hence, a mall value of GW can be 
considered a genetic signature of a population bottleneck 
(Excoffier and Lischer, 2010; Garza and Williamson, 
2002). The Garza-Williamson statistics were computed in 
the software Arlequin v3.5 (Excoffier and Lischer, 2010). 

Second, the signature of a population bottleneck was 
investigated based on heterozygosity excess (Cornuet 
and Luikartt, 1996), using Bottleneck v1.2.02 software 
(Luikart and Cornuet, 1999). In Bottleneck v1.2.02, 
three different mutation models, namely the infinite 
allele model (IAM), the stepwise mutational model 
(SMM), and the two-phase mutational model (TPM); and 
three different statistical tests namely the Sign test, the 
Standardized difference test and the Wilcoxon signed- 
rank test are available for heterozygosity excess testing 
(Luikart and Cornuet, 1999). Considering the number 
of used loci, heterozygosity excess was tested by the 
Wilcoxon signed-rank test, under all mutation models. 
The following settings were applied: TPM with SMM 
accounting for 95% of mutations, 5% SMM, a variance 
among multiple steps of 12 000 and 10 000 iterations 
(Luikart and Cornuet, 1999). 

Third, genetic signature of population bottleneck 
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was investigated by Mode-shift test (Luikart et al., 
1998), using the software Bottleneck v1.2.02 (Luikart 
and Cornuet, 1999). A putatively stable population is 
expected to have a peak of allele number at the lowest 
frequency class, resulting in an L-shaped distribution 
of allele frequency. Population bottleneck displaces the 
peak at other frequency classes, causing distortion of this 
graphical pattern or a shifted mode of allele frequency 
distribution (Luikart et al., 1998). 

Finally, the trajectory of the effective population size 
on long-term was investigated by comparing the historical 
and the contemporary values. Assuming mutation-drift 
equilibrium, the historical effective population size N, was 
estimated based on heterozygosity data as N, = (1/[1-H.]’- 
1)/8u, where H, is the average expected heterozygosity 
and и the mutation rate per generation (Nei, 1987). 
For и, 10°, the value reported in other amphibians was 
used (Peery et al., 2012). The contemporary N, was 
estimated based on LD, using the single-sample method 
as implemented in NeEstimator v2.0 software (Do et al., 
2014). A potential bias from rare alleles was controlled 
and a pcritic value of 0.01 was used (Do et al., 2014). 
LD estimates are based on the expectation that as №, 
decreases, genetic drift systematically associates alleles 
at different loci. Hence, in bottlenecked populations, the 
level of LD reflects №, (Do et al., 2014). 


2.4. Ethic statement Access to the focal species and 
the fieldwork did not require any formal research permit. 
The animal manipulation and sampling protocols were 
previously approved by the Animal Welfare Committee of 
the Central China Normal University. 


3. Results 


3.1. Adult census population size We found H. 
chinensis at 4 sites of different localities (Figure 1; Table 
1), in a total area of 2.18 кт”. These sites were located 
north-east of Yichang city (Figure 1). The overall adult 
census population size in the area slightly varied between 
breeding seasons. It ranged from 386 to 404 individuals 
(Table 1). Reproduction was completed in permanent 
ponds of varying sizes (Table 1). These ponds occurred 
at altitudes ranging from 1 408 to 1 582 m (Table 1) and 
were distant from 1.01 to 3.2 km. 

In Dongtoutang during the 2015-2016 breeding 
season, we identified 62 H. chinensis adult females and 
78 adult males, resulting in an N, of 140 individuals and 
a sex ratio of 1.26:1. We sampled 118 adult individuals in 
total, including 40 females and the 78 males. H. chinensis 
individuals appeared in the breeding pond found in this 
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locality on November 15. The number of male individuals 
found in the pond increased from the first (November 26, 
2015) to the third (January 4, 2016) census, reaching a 
peak of 60 individuals (Figure 2). The peak lasted till the 
sixth (February 20, 2016) census with minor fluctuation 
(Figure 2). The number of male individuals found in the 
pond then decreased till the last (March 8, 2016) census 
(Figure 2). 

As expected, we obtained a high recapture rate of male 
individuals (Figure 2). Out of the 78 male individuals 
identified in total, 67 (85.90%) were already identified 
and toe clipped on January 4, 2016, i.e. with the third 
census (Figure 2). Since the fourth census on January 
13, 2016, 97.30-10096 of all male individuals found 
in the breeding pond were individuals that had already 
been identified and toe clipped (Figure 2). As expected, 
together, this high recapture rate and the evolution pattern 
of the number of male individuals found in the pond 
along the breeding season suggest that male individuals 
spent several days in the breeding pond. They confirmed 
the accuracy of the method used to estimate the number 
of adult males and N.. 

The average value (+ SE) of N, over the study period 
(i.e. 2015-2018) varied between localities or breeding 
ponds. It ranged from 44 + 6 to 141 + 8 individuals (Table 
1). And, as expected, the maximum value of N, was 


observed in Dongtoutang (Table 1). The sex ratio was 
1.26: 1. 


3.2. Microsatellite markers We obtained 27 721 722 
(91.33%) clean reads in total. The Trinity assembly 
resulted in 77 111 unigenes. The size of the unigenes 
ranged from 201 to more than 3 000 bp (Figure 3). The 
total length of the unigenes and N50 were 46 797 297 bp 
and 951 bp, respectively. Among the unigenes, 46 378 
(61.1496) were not more than 400 bp in length. A total of 
20 599 (26.71%) unigenes were in the range of 401—1000 
bp. Unigenes longer than 3 000 bp were 1 507 (2.04%) in 
total (Figure 3). 

The transcriptomic sequencing raw data have been 
deposited into the NCBI Sequence Read Archive (SRA) 
database under the accession number SRP150396. They 
are accessible from https://www.ncbi.nlm.nih.gov/sra/ 
SRP150396. 

Microsatellite loci were identified in 3 516 unigenes. 
Among these unigenes, 2 999 (85.10%) contained 1 or 2 
microsatellite loci. Only 524 (14.90%) contained over 2 
microsatellite loci. We identified 4 297 microsatellite loci 
in total. The loci had repeat motifs varying from di- to 
tetra-. Among them, 3 254 (75.73%) loci had dinucleotide 
repeat motifs. Trinucleotide repeat motifs were observed 
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Figure 1 Locations of study sites of Hynobius chinensis in Changyang county. Sites are distinguished by the names of the corresponding 


localities. DNG: Dongtoutang; SHJ: Shuijingwan; JNJ: Jiangjunao; HNG: Hengtun. 
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Figure 2 Number of Hynobius chinensis male individuals present in a breeding pond at different dates in 2015—2016 in Dongtoutang 


(Changyang county, Hubei). 


Table 1 Adult census population size of Hynobius chinensis estimated at different breeding ponds in different localities in Changyang 


County. 

Locality Longitude Latitude Elevation (m) 
Dongtoutang 110°41'49" E 30?34'42" N 1548 
Shuijingwan 110°41'24" E 30?35'04" N 1582 
Jiangjunao 110°40'25" E 30°34'59" N 1482 
Hengtun 110°39'34" E 30?35'16" N 1408 


PS (m) N, 
7 14148 
25 91412 
20 44+6 
15 111 4 10 


№: adult census population size; PS: pond size. For each sex, and for N,, values are given in terms of average (+ standard error) over the 


study period, i.e. 2015-2018. 
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in 223 (5.2%) loci. Tetranucleotide repeat motifs were 
observed in 120 (2.7994) loci. 

We developed and characterized 13 polymorphic 
microsatellite loci. These loci were in HWE and there 
were no evidence of significant LD between pairs of 
loci, of null alleles or genotyping errors in the 22 sample 
dataset. Among the 13 microsatellite loci, 7 (54%) loci 
have tetranucleotide repeat motifs and 6 (46.15%) have 
dinucleotide repeat types. The remaining locus has a 
trinucleotide repeat motif. The number of repeats range 
from 5 to 9 and the length of the loci range from 105 to 
253 bp (Table 2). 

The microsatellite DNA sequences have been deposited 
in GenBank. The characteristics and the accession number 
of each locus are shown in Table 2. 


3.3. Genetic diversity and signature of a population 
bottleneck Genotype data at the 13 microsatellite loci 
have been obtained for all the 118 sampled H. chinensis 
adult individuals (see Table S1 in Appendix). As observed 
during the development of the microsatellite markers, the 
analyses in Micro-checker v2.2.3 provided no evidence 
of either null alleles or genotyping errors. And, there 
were no cases of linkage disequilibrium among the loci. 
Among the 13 loci, the locus HC64 revealed a relatively 
low PIC (Table 2) and 5 loci, i.e. HC23, HC47, HC54, 
HC85 and HC116 significantly deviated from HWE, after 
Bonferroni correction (P < 0.0042). All these 5 loci and 
the locus HC64 were therefore removed from further 
analyses. As a result, 7 microsatellite loci (Table 3) were 
used in the analyses. 

For the 118 samples and 7 loci, 15 different alleles were 
identified, resulting in 2.14 alleles per loci or average A 
of 2.14 (Table 3). In other terms, the value of A was 2 for 
6 (85.71%) out of the 7 loci (Table 3). Allele frequency 
ranged from 0.013 to 0.87. The maximum number of 
alleles (46.15%) was observed for the frequency range of 
0.4—0.6 (Figure 4). 23.07% of alleles were found in the 
lowest frequency range, 1.e. 0—0.2 (Figure 4). Average 
H, and H, were 0.38 and 0.40, respectively. The average 
value of the inbreeding coefficient, F, was —0.06 (Table 
3). 

The value of GW varied between loci, ranging from 
0.17 to 1, and with an average value of 0.61 (Table 3). 
Comparison of each locus’GW to the critical value for 
rejecting the null hypothesis of demographic stability 
(i.e. 0.68) showed that GW was below this value for 4 
(57%) loci (Table 3). The Wilconxon signed-rank test 
of heterozygosity excess yielded significant (One tail 
for heterozygosity excess) results (P « 0.05) as well as 
under IAM (P = 0.004), SMM (P = 0.004) and STM 
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Figure 3 Distribution of unigenes by size range in the transcriptome 
of Hynobius chinensis. 
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Figure 4 Distribution of alleles at seven microsatellite loci by 
frequency range in Hynobius chinensis (N — 118). 


(P — 0.027). In the Sing test, heterozygosity excess was 
detected for all the 7 loci under IAM. Under TPM and 
also under SMM, it was detected for 5 loci namely HC19, 
HC26, HC33, HC103 and HC114. A significant distortion 
of the typical L-shaped distribution of allele frequency of 
putatively stable populations was detected. The historical 
N, was estimated at 234 individuals. The contemporary N, 
was found to be 36 (95% CI: 13.1—124.6) individuals. 


4. Discussion 


With an average value of 0.61, the results of Garza- 
Williamson test indicate that the studied population of H. 
chinensis has experienced a severe reduction in size or 
population bottleneck. The results of Wilconxon signed- 
rank test and Mode shift test also evidence signatures of a 
population bottleneck. 

Detection of a recent population bottleneck in H. 
chinensis seems to not a surprising finding (Fei and 
Ye, 2016; Wang et al., 2007). Literature shows that H. 
chinensis may have experienced a severe demographic 
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Table3 Values of genetic diversity indices and Garza-Williamson statistics (N — 118) 


Locus А ТЇ; А F,, Size range (bp) GW 
НС19 2 0.318 0.384 0.143 182—198 0.4 
НС28 2 0.591 0.511 0.042 133-177 0.167 
HC33 2 0.591 0.485 —0.001 173-177 1 
НС49 2 0.136 0.206 0.102 249—253 1 
HC63 3 0.273 0.312 —0.152 105-111 1 
HC103 2 0.364 0.304 —0.147 212—222 0.333 
НС114 2 0.409 0.511 0.009 161—169 0.4 
Mean 2.14 0.383 0.395 —0.06 0.614 


A: number of alleles; H,: observed heterozygosity; H,: expected heterozygosity; F;,: inbreeding coefficient; GW: Garza-Williamson index. 


crisis in the last two centuries (Wang et al., 2007). 
For example, Wang et al. (2007) reported that after its 
description at the end of the 19" century, H. chinensis 
became irretrievable till 2005 in the current location of 
the Yichang city (Wang et al., 2007). 

With no more than 2.14 alleles per locus, and with 
0.38 and 0.40 as respective average values of H, and 
H,, the genetic diversity revealed by H. chinensis 
seems obviously low in comparison to that reported in 
Least Concern (IUCN, 2018) pond-breeding Hynobiid 
salamanders (Matsunami et al., 2015; Yoshikawa et al., 
2013). For examples, the study of small samples from 
2 populations (№ = 16 and N = 8) revealed 6.63 and 
3.53 alleles per locus in Н. nebulosus (Yoshikawa et al., 
2013). Another study of small samples (N — 20) from 
two populations of H. retardatus found that the average 
numbers of alleles per locus were 4.25 and 3.83. The 
average H, for these populations was 0.64 and 0.57; and 
the average H, was found to be 0.59 and 0.49 (Matsunami 
et al., 2015). 

The level of genetic diversity revealed by H. chinensis 
in the present study seems close to that reported in other 
threatened and narrowly distributed (IUCN, 2018) closely 
related species such as H. dunni (Sugawara et al., 2015), 
H. maoershanensis (Lin et al., 2015), and H. tokyoensis 
(Sugawara et al., 2015). A recent study of 12 populations 
of Н. dunni reported an average А of 2.53 and 0.34 as 
the value of H, and H, (Sugawara et al., 2015). In a 
another recent large-scale study of H. tokyoensis, a forest- 
dwelling Hynobiid salamander threatened by habitat loss 
and fragmentation due to land use changes as it might be 
in the focal species (Fei and Ye, 2016) which involved 
46 populations from 12 different geographic regions, 
Sugawara et al. (2016) found comparable results. These 
authors showed that the maximum average value of А 
was 3.4. They also showed that average А was « 2.0 
in 39.13% of populations and between 2.0 and 2.40 in 
32.61% of populations (Sugawara et al., 2016). H, ranged 
from 0.02 to .51 and averaged 0.27 for all the populations. 


Our demographic and genetic diversity results seem to 
corroborate findings of earlier studies, using microsatellite 
data, which report population bottlenecks (Apodaca et al., 
2012; Heredia-Bobadilla et al., 2017; Wang et al., 2017), 
and levels of genetic diversity similar to that revealed by 
Н. chinensis in this study (Jordan et al., 2009; Spear et al., 
2006) in a variety of threatened salamanders, including 
Ambystomatids, and Phlethodontids. 

Inbreeding is one of the mechanisms expected to lead 
to loss of genetic diversity after population bottlenecks in 
threatened and narrowly distributed amphibians (Ficetola 
et al., 2011; Storfer et al., 2014). However, with an 
average value of —0.06, the inbreeding coefficient, F; 
provides evidence of moderate level of inbreeding in the 
present study. 

Analyses of microsatellite data in Genepop 
v1.2 software detected departures from HWE for 5 
microsatellite loci (i.e. 2C23, HC47, HC54, HC85 and 
НС116). These departures cannot be attributed to null 
alleles because a significant occurrence of such alleles in 
the dataset has not been confirmed by analyses in Micro- 
Checker v2.2.3. In comparison to loci in HWE, these loci 
exhibited heterozygosity deficiency (Results not shown). 
Hence, though inbreeding coefficient F, estimated without 
these loci revealed a moderate level of inbreeding, the 
departures from HWE can be considered as a consequence 
of inbreeding following the inferred population bottleneck 
(Apodaca et al., 2012; Storfer et al., 2014). Inbreeding 
may have caused loss of heterozygosity, which resulted 
in the observed low level of heterozygosity (Storfer et al., 
2014). The level of inbreeding may have been brought at 
a moderate level by selection against inbreeding (Ficetola 
et al., 2011) and by polyandry (Slatyer et al., 2012). Few 
inbreed individuals may reach maturity (Ficetola et al., 
2011). During this study, at all observation opportunities 
(8 in total), we found that an individual female’s eggs 
were seemingly fertilized by several males (2-7). Such 
a mating pattern may result in a reduced probability of 
inbreeding (Slatyer et al., 2012). 
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The low number of alleles revealed by H. chinensis in 
this study may have been a consequence of loss of rare 
alleles by genetic drift following the inferred population 
bottleneck (Garza and Williamson, 2001). Loss of such 
alleles may have been facilitated by fragmentation of 
populations caused by forced exclusion of H. chinensis 
from proportions of natural habitats (Allentoft and 
O’Brien, 2010; Apodaca et al., 2012; Sugawara et al., 
2016). Also, a low density of breeders (i.e. 386—404 
breeders per 2.18 km’) may have contributed to 
fragmentation of populations (Allentoft and O’Brien, 
2010). 

According to interviewed people, environmental 
conditions have considerably changed in the study area 
since 1960s. Since that time period, the local forest that is 
also the habitat of H. chinensis was intensively harvested. 
From valleys toward mountain tops, cultivated lands 
expanded rapidly, reaching almost the current area in 
1970s. The forest rapidly regressed towards mountain tops 
and became fragmented. As larger trees were selectively 
cut, the density of the canopy rapidly decreased and the 
forest become degraded. Cultivated lands rapidly became 
infertile, and the use of fertilizers became intensive in 
1970s. 

Interviewed people indicated that fire was largely used 
to clear forest areas for land cultivation, suggesting a 
possible massive and repeated decimation of terrestrial 
life stages of H. chinensis. Human-induced forest 
reduction, fragmentation and degradation were evident 
in all the study localities. Exposure of H. chinensis (both 
aquatic and terrestrial stages) to agrochemicals was also 
suspected. 

These human-induced changes in the potential habitat 
may have affected gene flow and population growth 
rate and caused population fragmentation and declines, 
triggering loss of alleles and heterozygosity in H. 
chinensis (Allentoft and O’Brien, 2010; Apodaca et al., 
2012; Heredia-Bobadilla et al., 2017; Rivera-Ortíz et al., 
2015; Sugawara et al., 2016). 

In addition, interviewed people revealed that during 
breeding seasons, adult H. chinensis individuals are 
collected from breeding ponds and used as a source of 
proteins. This finding suggests that overharvesting may 
be another contributing factor to population declines and 
loss of genetic diversity (Hauser et al., 2002) in the study 
system. 

The values of adult census population size, N, revealed 
by H. chinensis in the present study appears to be smaller 
than values reported by earlier studies of close related 
species under suitable habitat conditions (Fu et al., 2003; 


Gu et al., 1999; Ma and Gu, 1999), suggesting a possible 
reduced abundance of H. chinensis in the study area due 
to threats mentioned above (Chen et al., 2016; Fei and 
Ye, 2016). For example, survey data from H. amjiensis 
gathered from 1992 to 1998 (Gu et al., 1999), from H. 
yiwuensis -identified by the authors as H. chinensis- 
gathered in 1985, 1988 and 1998 (Ma and Gu, 1999) and 
personal observation of H. yiwuensis by Fu et al. (2003) 
show that in natural and less anthropogenically disturbed 
habitats, the number of female breeders per breeding 
pond should exceed 50, even 100 individuals (Chen et al., 
2016; Fu et al., 2003; Gu et al., 1999; Ma and Gu, 1999). 

Comparison of the contemporary N, to the historical 
М, shows that the contemporary N, is more than 6 orders 
of magnitude smaller than the historical N,, evidencing a 
long-term reduction in population size. This suggests that 
H. chinensis may have undergone population declines and 
loss of genetic diversity prior to the anthropogenic threats 
listed above (Jordan et al., 2009). And, long-term declines 
may have also significantly contributed to fragmentation 
in gene pools and loss of genetic diversity (Jordan et al., 
2009). 

We report data on М, in Н. chinensis from three 
different localities collected over three successive 
breeding period. Despite a relative stability during the 
study period, as №, may naturally highly fluctuate, long- 
term data may be required to know the exact №, in the 
study localities (Beebee and Griffiths, 2005; McCartney- 
Melstad and Shaffer, 2015). The observed variability in N, 
between localities suggests that local factors may govern 
the distribution and abundance of H. chinensis (Beebee 
and Griffiths, 2005; Gu et al., 1999). And, this should 
also be true for genetic diversity (Beebee and Griffiths, 
2005; Sugawara et al., 2016). Evidence from closely 
related species such as H. amjiensis shows that variation 
in №. may be caused by spatial variation in human- 
induced changes in ponds’ water quality, and physical 
characteristics (Gu et al., 1999). 

We sampled the largest population found, i.e. 
Taungtotang's population expecting to see a better genetic 
variability and demographic history (Hoffmann et al., 
2017). But our results show that despite its bigger size, 
this population is genetically poor and bottlenecked. 
And, as expected for amphibians (Allentoft and O'Brien, 
2010; Beebee and Griffiths, 2005), N, was several 
orders of magnitudes smaller than N,. Note that, in 
2015, we recorded 140 H. chinensis adult individuals 
in Taungtotang, resulting in a N/N, ratio of 0.26 (or 
36:140). Though literature shows that higher values can 
be reached in some species (e.g. Georcrina vitellina, 
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Rana sylvatica), the value of this ration is « 0.2 in many 
amphibians, including anurans and Caudata (Beebee and 
Griffiths, 2005). As N, was smaller in the other localities 
than in Taungtotang, this relationship suggests that N, and 
genetic diversity must be more reduced in those localities 
(Hoffmann et al., 2017), suggesting a possible extirpation 
of H. chinensis’ populations. 

As mentioned above, information on landscape 
use lacks for H. chinensis. The results of this study 
demonstrate the paucity of the existing information. For 
example, while expected to occur at altitudes ranging 
from 1400 to 1500 m (Fei and Ye, 2016), we found 
two (5096) breeding populations, including the largest 
breeding populations beyond this range, at altitudes 
of 1548 m and 1582 m precisely. The lack of such 
information prevents us from appreciating the sampling 
effort of this study and hence the significance of the 
results at the species level (Allentoft and O'Brien, 2010; 
Beebee and Griffiths, 2005). 

In conclusion, this study makes available thirteen 
novel microsatellite markers for H. chinensis. It provides 
information on the exact location of populations and 
breeding sites. This study provides molecular evidence of 
a population bottleneck in H. chinensis, which is provided 
by the Garza-Williamson test, Heterozygosity excess 
tests, Mode shift test and estimates of effective population 
sizes. It shows that H. chinensis may exist in small-sized 
populations, which retain low microsatellite genetic 
diversity. Together, our results suggest that populations of 
H. chinensis may have been extirpated in the study area. 

For further understanding of the demography and 
genetic viability of H. chinensis, future studies should also 
investigate the distribution and sizes of the representative 
populations of this endangered species. They should 
investigate fitness effects of inbreeding and specific 
causes of population declines while considering gene 
flow and populations genetic structuring. Conservation 
initiatives should aim to conserve and restore the habitats 
and regulating the exploitation of the species. They 
should aim to maintain the extant genetic variation. 
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